Module 1 Project - King County House Sales¶
Matt Sparr - Self-Paced Data Science Program - Review 10/31/2018¶
Introduction¶
The purpose of this project is to investigate the variables of homes sold in King County, Washington and to develop a model to accurately predict home sale prices. We will be following the OSEMiN model which involves the steps of obtain, scrub, explore, model, and finally interpret.
Obtaining the Data¶
We will first import all the necessary libaries and modules and then use Pandas to grab the dataset from the kc_house_data.csv and store it as a Pandas Dataframe.
In [73]:
#importing required libraries and setting matplotlib to inline
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import math
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression
import sklearn.linear_model as lm
from sklearn.feature_selection import RFE
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from statsmodels.formula.api import ols
%matplotlib inline
In [2]:
#grabbing the dataset
df = pd.read_csv('kc_house_data.csv')
Scrubbing the Data¶
Now we will clean the data to enable us to create our model without errors. First we will drop the 'id', 'date', 'lat', 'long', and 'zipcode' columns as they will not be needed. We will then inspect the dataframe and address the columns will null or missing values.
In [3]:
#dropping id and date columns
df.drop(['id', 'date', 'lat', 'long', 'zipcode'], axis=1, inplace=True)
In [4]:
#viewing the first 15 rows
df.head(15)
Out[4]:
In [5]:
#checking value in sqft_basement
df.sqft_basement.value_counts()
Out[5]:
There appears to be some values of '?' in sqft_basement. We will replace these with the column mean.
In [6]:
#replacing missing data in sqft_basement with column mean
sqft_basement_numerical = df.sqft_basement[df.sqft_basement != '?']
sqft_basement_numerical = sqft_basement_numerical.astype(float)
sqft_basement_numerical_mean = round(sqft_basement_numerical.mean(),1)
df['sqft_basement'] = (df['sqft_basement'].map(lambda x: sqft_basement_numerical_mean if x == '?' else x)).astype(float)
In [7]:
#checking for null values
df.isna().sum()
Out[7]:
The 'view', 'waterfront', and 'yr_renovated' columns all have null values. We will assume that those homes had no views, are not waterfront, and/or were not renovated and fill those null values with 0.
In [8]:
#replacing null values in waterfront with 0
df['view'] = df['view'].fillna(0)
df['waterfront'] = df['waterfront'].fillna(0)
df['yr_renovated'] = df['yr_renovated'].fillna(0)
Now before we modify any of the columns we will create a copy of the dataframe to give us an unaltered version to run models on later for comparison.
In [9]:
#creating copy of df before transformations
df_original = df.copy()
Exploring the Data¶
For this section, we will begin by inspecting the histogram of each variable to look at the distributions and skew.
In [10]:
#intial histogram of all variables
df.hist(figsize=(15,15))
Out[10]:
An initial observation of the histograms reveals that nearly none of the variables exhibit a normal distribution and there are outliers in many of the variables. To address these issues, we will iterate through each variable and perform adjustments and transformations in order to give us a more normal distribution.
In [11]:
#checking bedrooms for outliers
df.bedrooms.value_counts()
Out[11]:
In [12]:
#histogram before removal of outliers
df.bedrooms.hist()
Out[12]:
In [13]:
#removing outliers
df.drop(df.index[df['bedrooms'] >= 8], inplace=True)
In [14]:
#histogram after removal of outliers
df.bedrooms.hist()
Out[14]:
The 'bedrooms' column had some outliers so all homes with more than 7 bedrooms were removed from the dataframe. The number of homes removed was 24 which represents an insignificant percentage of our total data so removing them should not significantly impact the final model.
In [15]:
#bathrooms histogram
df.bathrooms.hist()
Out[15]:
In [16]:
#using log transformation on bathrooms
(np.log(df.bathrooms)).hist()
Out[16]:
In [17]:
#transforming bathrooms
df.bathrooms = np.log(df.bathrooms)
The 'bathrooms' column showed a negatively skewed distribution so we transformed it using a log transformation as the data now exhibits a more normal distribution.
In [18]:
#sqft_living histogram
df.sqft_living.hist()
Out[18]:
In [19]:
#checking count of outliers
print(len(df[df.sqft_living > 8000]))
print(len(df[df.sqft_living > 6000]))
print(len(df[df.sqft_living > 4000]))
In [20]:
#removing outliers
df.drop(df.index[df['sqft_living'] > 4000], inplace=True)
In [21]:
#histogram after removal of outliers
df.sqft_living.hist()
Out[21]:
In [22]:
#histogram of square root transformed sqft_living
(np.sqrt(df.sqft_living)).hist()
Out[22]:
In [23]:
#square root transformation on sqft_living
df.sqft_living = np.sqrt(df.sqft_living)
The 'sqft_living' column had some significant outliers. Homes with more than 4000 square feet were removed which was a little over 800 entries. The histogram of the modified data still showed some negative skew so we then performed a square root transformation thus giving us a more normal distribution.
In [24]:
#sqft_above histogram
df.sqft_above.hist()
Out[24]:
In [25]:
#square root transformation
df.sqft_above = np.sqrt(df.sqft_above)
In [26]:
#square root transformed histogram
df.sqft_above.hist()
Out[26]:
The 'sqft_above' column showed a negative skew so we transformed it using a square root transformation to make the distribution more normal.
In [27]:
#sqft_basement histogram
df.sqft_basement.hist()
Out[27]:
In [28]:
#creating new column that declares whether or not a home has a basement
df['basement'] = df.apply(lambda row: 0 if row.sqft_basement == 0 else 1, axis=1)
In [29]:
#new binary basement histogram
df.basement.hist()
Out[29]:
The 'sqft_basement' column has a lot of values of 0 which heavily skews the distribution. By creating a new column 'basement' which is binary in nature and declares whether or not a home has a basement, we will better be able to see how this variable affects the home price later on in our model.
In [30]:
#viewing stats of yr_built
df.yr_built.describe()
Out[30]:
In [31]:
#creating new column age that declares the age of the home
df['age'] = (df.yr_built.max() - df.yr_built + 1)
In [32]:
#histogram of new age column
df.age.hist()
Out[32]:
In [33]:
#square root transformation on age
df.age = np.sqrt(df.age)
df.age.hist()
Out[33]:
In [34]:
#removal of old yr_built column
df.drop('yr_built',axis=1,inplace=True)
The 'yr_built' column did not make a lot of sense to have in its current state as it represents the age of the home but in terms of the year built. To remedy this we created a new column 'age' which subtracts the year the home was built from the most recent year in the data - 2015. A value of 1 was also added to the age in order to allow us to transform it without error.
In [35]:
#sqft_lot histogram
df.sqft_lot.hist()
Out[35]:
In [36]:
#checking number of outliers
len(df[df.sqft_lot>200000])
Out[36]:
In [37]:
#removal of outliers
df.drop(df.index[df['sqft_lot'] > 200000], inplace=True)
In [38]:
#histogram after removal of outliers
df.sqft_lot.hist()
Out[38]:
In [39]:
#log transformation of sqft_lot
df.sqft_lot = np.log(df.sqft_lot)
df.sqft_lot.hist()
Out[39]:
For the 'sqft_lot' column, we first had to remove some of the outliers (213 values) and then performed a log transformation on the data to normalize the distribution.
In [40]:
#histogram of sqft_lot15
df.sqft_lot15.hist()
Out[40]:
In [41]:
#log transformed sqft_lot15
df.sqft_lot15 = np.log(df.sqft_lot15)
df.sqft_lot15.hist()
Out[41]:
Similarly to the 'sqft_lot' column, we performed a log transformation on the 'sqft_lot15' column to normalize the distribution.
In [42]:
#view histogram
df.view.hist()
Out[42]:
In [43]:
#unique values of view
df.view.unique()
Out[43]:
In [44]:
#converting all view values above 0 to 1 to signify that the home has been viewed
df['viewed'] = df.view
i = 4
while i > 0:
df.viewed = df.viewed.replace(i, 1)
i -= 1
df.viewed.hist()
Out[44]:
The 'view' column only contains values 0, 1, 2, 3 and 4 with most of the values being 0. To better represent this data, an additional column was added which tells us whether or not a home was view, represented by a 1 or 0.
In [45]:
#price histogram
df.price.hist()
Out[45]:
In [46]:
#log transformation on price
df['price_log'] = np.log(df.price)
df.price_log.hist()
Out[46]:
Finally, we perform a log transformation on our dependent variable 'price' in order to normalize its distribution.
In [47]:
#updated histogram of all variables after transformations and adjustments
df.hist(figsize=(14,14))
Out[47]:
Above we can see all of the updated histograms for our variables. Many of the distributions are now normal apart from our created binary columns and their original counterparts.
Model the Data¶
Now it is time to create a model for our data. We will be using the function 'stepwise_selection' below to rank the predictors of our data on their contribution to the model and to find the best number of predictors to use.
In [48]:
#stepwise selection function
import statsmodels.api as sm
def stepwise_selection(X, y,
initial_list=[],
threshold_in=0.01,
threshold_out = 0.05,
verbose=True):
""" Perform a forward-backward feature selection
based on p-value from statsmodels.api.OLS
Arguments:
X - pandas.DataFrame with candidate features
y - list-like with the target
initial_list - list of features to start with (column names of X)
threshold_in - include a feature if its p-value < threshold_in
threshold_out - exclude a feature if its p-value > threshold_out
verbose - whether to print the sequence of inclusions and exclusions
Returns: list of selected features
Always set threshold_in < threshold_out to avoid infinite looping.
See https://en.wikipedia.org/wiki/Stepwise_regression for the details
"""
included = list(initial_list)
while True:
changed=False
# forward step
excluded = list(set(X.columns)-set(included))
new_pval = pd.Series(index=excluded)
for new_column in excluded:
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included+[new_column]]))).fit()
new_pval[new_column] = model.pvalues[new_column]
best_pval = new_pval.min()
if best_pval < threshold_in:
best_feature = new_pval.idxmin()
included.append(best_feature)
changed=True
if verbose:
print('Add {:30} with p-value {:.6}'.format(best_feature, best_pval))
# backward step
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
# use all coefs except intercept
pvalues = model.pvalues.iloc[1:]
worst_pval = pvalues.max() # null if pvalues is empty
if worst_pval > threshold_out:
changed=True
worst_feature = pvalues.argmax()
included.remove(worst_feature)
if verbose:
print('Drop {:30} with p-value {:.6}'.format(worst_feature, worst_pval))
if not changed:
break
return included
We will now split our data into a train and test set using an initial test size of 0.25. We will then use the stepwise function above to do an interative feature selection for all of the independent variables. For each iteration, one more predictor will be selected and added, and in addition, a repeated k-fold with 3 splits and 10 repeats will be ran to cross validate the R-squared value and % difference in MSE between the train and test set for each regression. Then, by graphing a plot of our predicted y-values vs the actual y-values and comparing the R-squared values and % differences in MSE, we should be able to decide which predictors to include in our model of best fit.
In [49]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df.drop(['price', 'price_log'], axis=1)
y = df.price
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df.drop(['price', 'price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
#creating a plot of predicted values vs actual values
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df['# Predictors'] = predictors
results_df['R-squared'] = reg_score
results_df['MSE % Difference'] = mse_diffs
results_df['Added Predictor'] = added_pred
plt.show()
results_df
Out[49]:
From the above table, we can see that the test set MSE is anywhere from approx 1% more to 1% less than the training set MSE. This varies slightly as more predictors are added and doesn't follow a pattern unlike the R-squared values which increases with each additional predictor. Looking at the generated plots of measured vs predicted values for 'price', we can see that as more predictors are added, the predicted values approach the measured values. Based on these findings, we will use all 17 predictors in our model as that will maximize our R-squared value of the model against our test data and also will maintain a moderately low difference in MSE between the train and test set.
Evaluate the Model¶
Now that we have our model, we can begin to ask relevant questions as to the effectiveness and fit of the model as well as address concerns and investigate possible improvements.
Before we started to dig into and modify our data set, we grabbed a copy of an unedited version of the data - 'df_orig'. We can run a regression on that dataset and compare the results to those of our model to answer this question. We will repeat the stepwise feature selection as we did earlier, using a repeated k-fold validation method, and plotting the predicted values vs the actual values.
In [50]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_original.drop(['price'], axis=1)
y = df_original.price
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_original = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_original.drop('price', axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_original['# Predictors'] = predictors
results_df_original['R-squared'] = reg_score
results_df_original['MSE % Difference'] = mse_diffs
results_df_original['Added Predictor'] = added_pred
plt.show()
results_df_original
Out[50]:
Comparing the regression of the data in 'df_original' to that of our model we can see that similarly, the R-squared value of the test set increases as the number of predictors increase apart from the last step in going from 13 to 14 predictors where the R-squared value slightly decreases. The R-squared of our model using 17 predictors was 0.614463 while the best R-squared of the un-transformed dataset is actually higher at 0.652061 using 13 predictors. This is interesting as based on that metric alone, the model of 'df_original' actually fits the test set slightly better on average. Conversely, however, the % difference in MSE between the train and test set is lower in our model with it being -0.943803 and for 'df_original' using 13 predictors, 5.093368.
When we look at the plots of predicted vs measured values, we can see that the model of 'df_original' fits the data moderately well as more predictors are added. However, the model of 'df_original' starts to do a worse job of predicting as the 'price' increases.
It would appear that our model does fit the data better and is a better predictor of home prices, but to a lesser extent than would be expected. To investigate why this might be the case, we can ask further questions.
If multiple variables exhibit covariance with each other, our model could be affected negatively. We will check the covariance of all of our predictors and try eliminating problematic variables to see if it improves our model.
In [51]:
#heatmap to check for correlations between variables where corr > 0.75
sns.heatmap(df.corr(),center=0);
Some of our variables appear to be highly correlated. We will choose a threshold of 0.75 for variables with significant correlation and generate a new heatmap to more clearly show those variables.
In [52]:
sns.heatmap(abs(df.corr())>0.75, center=0)
Out[52]:
From the above we can see that 'sqft_above' and 'sqft_living', 'sqft_lot15' and 'sqft_lot', 'view' and 'viewed', and 'sqft_basement' and 'basement', show high correlation with each other. Out of each of these pairs, the less important ones are 'sqft_above' as the total square footage of a home is likely more significant, 'sqft_lot', 'viewed', and 'basement'. Now let's try running our regression after removing these variables.
In [53]:
#creating dataframe that will remove some of the variables with high correlation to other variables
df_adj = df.copy()
to_drop = ['view', 'sqft_lot15', 'sqft_above', 'sqft_basement']
df_adj.drop(to_drop, axis=1, inplace=True)
In [54]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_adj.drop(['price','price_log'], axis=1)
y = df_adj.price_log
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_adj = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_adj.drop(['price','price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_adj['# Predictors'] = predictors
results_df_adj['R-squared'] = reg_score
results_df_adj['MSE % Difference'] = mse_diffs
results_df_adj['Added Predictor'] = added_pred
plt.show()
results_df_adj
Out[54]:
Looking at the results of this regression, we can see that the R-squared value and % difference in MSE are approximately the same in this adjusted model as in our original model. This adjusted model is a slightly worse fit as the R-squared values are slightly lower and the % differences in MSE are slightly more. This seems to suggest that the covariance of the removed variables does not significantly affect our model's fit in a negative way. One additional thing to consider however, is our added variables of basement and viewed. Both of which are binary data columns meant to deal with both columns having a large number of 0s and then a small spread of other values. We will try running a regression simply without those two added variables to see if that makes any difference.
In [55]:
#creating dataframe that will remove some of the variables with high correlation to other variables
df_adj2 = df.copy()
to_drop = ['viewed', 'basement']
df_adj2.drop(to_drop, axis=1, inplace=True)
In [56]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_adj2.drop(['price', 'price_log'], axis=1)
y = df_adj2.price_log
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_adj2 = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_adj2.drop(['price', 'price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_adj2['# Predictors'] = predictors
results_df_adj2['R-squared'] = reg_score
results_df_adj2['MSE % Difference'] = mse_diffs
results_df_adj2['Added Predictor'] = added_pred
plt.show()
results_df_adj2
Out[56]:
Again, it seems that the fit is similar in goodness to our model but slightly worse as marked by slightly higher R-squared values and slightly more % difference in MSE between the test and train sets. For this reason we will keep our model as is, including our created columns of 'basement' and 'viewed'.
Let's begin by viewing the statistics of 'price' and its distribution.
In [57]:
df.price.describe()
Out[57]:
In [58]:
df.price.hist()
Out[58]:
75% of homes have a 'price' of 615,000 or less. Let's try using this as our subset and seeing how that affects our generated models.
In [59]:
df_s = df[df.price <= 615000]
In [60]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_s.drop(['price', 'price_log'], axis=1)
y = df_s.price_log
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_s = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_s.drop(['price', 'price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_s['# Predictors'] = predictors
results_df_s['R-squared'] = reg_score
results_df_s['MSE % Difference'] = mse_diffs
results_df_s['Added Predictor'] = added_pred
plt.show()
results_df_s
Out[60]:
Immediately we can see that this did not work. While the % difference in MSE is low, and lower than our original model, the R-squared value was nearly cut in half. Maybe a different subset will work. We will try both above and below the 50% mark of 440,000.
In [61]:
df_s2 = df[df.price <= 440000]
In [62]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_s2.drop(['price', 'price_log'], axis=1)
y = df_s2.price_log
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_s2 = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_s2.drop(['price', 'price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_s2['# Predictors'] = predictors
results_df_s2['R-squared'] = reg_score
results_df_s2['MSE % Difference'] = mse_diffs
results_df_s2['Added Predictor'] = added_pred
plt.show()
results_df_s2
Out[62]:
In [63]:
df_s3 = df[df.price >= 440000]
In [64]:
#RFE regression of all predictors using repeated k-fold
#setting X, y and creating train/test split
X = df_s3.drop(['price', 'price_log'], axis=1)
y = df_s3.price_log
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state=47)
#initializing counter variable 'i' and creating lists to add to over each iteration
i = 1
results_df_s3 = pd.DataFrame()
predictors = list()
reg_score = list()
mse_diffs = list()
added_pred = list()
previous_columns = []
#this loop repeats once for each column in 'df'
while i <= (int(len(df_s3.drop(['price', 'price_log'], axis=1).columns))):
rkf = RepeatedKFold(n_splits=3, n_repeats=10, random_state=47)
scores = 0
percent_diffs = 0
#this loop using a repeated k-fold to generate an average R-squared and % difference in train/test MSE
for train_index, test_index in rkf.split(X_train):
X_train2, X_val = X_train.iloc[train_index], X_train.iloc[test_index]
y_train2, y_val = y_train.iloc[train_index], y_train.iloc[test_index]
linreg = LinearRegression()
selector = RFE(linreg, n_features_to_select = i)
selector = selector.fit(X_train2, y_train2)
selected_columns = X_train2.columns[selector.support_ ]
linreg.fit(X_train2[selected_columns],y_train2)
scores = scores + linreg.score(X_test[selected_columns], y_test)
y_hat_train = linreg.predict(X_train[selected_columns])
y_hat_test = linreg.predict(X_test[selected_columns])
train_mse = mean_squared_error(y_train, y_hat_train)
test_mse = mean_squared_error(y_test, y_hat_test)
mse_diff = test_mse - train_mse
percent_diff = mse_diff/train_mse
percent_diffs += percent_diff
#adding data to lists
predictors.append(i)
reg_score.append(scores/30)
mse_diffs.append(percent_diffs/30*100)
added_pred.append(list(set(selected_columns) - set(previous_columns)))
previous_columns = selected_columns
predicted = linreg.predict(X_test[selected_columns])
fig, ax = plt.subplots(figsize=(5,5))
#ax.scatter(y_test, predicted, edgecolors=(0, 0, 0))
ax.plot(y_test, y_test, 'k--', lw=4, color='y')
sns.regplot(x=y_test, y=predicted, ax=ax, line_kws={"color": "red"})
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
ax.set_title(str(i) + " Predictors")
ax.legend()
i += 1
results_df_s3['# Predictors'] = predictors
results_df_s3['R-squared'] = reg_score
results_df_s3['MSE % Difference'] = mse_diffs
results_df_s3['Added Predictor'] = added_pred
plt.show()
results_df_s3
Out[64]:
Although choosing a subset where 'price' is greater than or equal to 615,000 gave the best result out of these experiments. All of them perform worse than our original model and signify that it is best we include the full range of home prices.
Interpreting the Data¶
In this final section we will use our model to draw conclusions about the data and its real-world applications.
We will look at the coefficients of our regression model. Since our dependent variable 'price' was log transformed, we will reverse that transformation for each of the coefficients in order to be able to talk about them in terms of real-life values. We will also sort the coefficients in descending order to rank the significance of each predictor.
In [65]:
#creates a dataframe of both log transformed and non-log transformed coefficients of predictors in model
X = df.drop(['price', 'price_log'], axis=1)
y = df.price_log
model_reg = LinearRegression()
linreg.fit(X, y)
model_coef = pd.DataFrame()
preds = list()
coef_log = list()
coef = list()
sum_coef = 0
percent_coef = list()
i = 0
for col in X.columns:
preds.append(X.columns.values[i])
coef_log.append(linreg.coef_[i])
coef.append(10**(linreg.coef_[i]))
sum_coef += 10**(linreg.coef_[i])
i += 1
model_coef['Predictor'] = preds
model_coef['Coefficient(log)'] = coef_log
model_coef['Coefficient'] = coef
for c in model_coef.Coefficient:
percent_coef.append(c/sum_coef*100)
model_coef['Contribution(%)'] = percent_coef
model_coef.sort_values(by=['Coefficient'], ascending=False)
Out[65]:
From the above table we can see the relative impact of each predictor with 'waterfront' being the most significant and 'sqft_lot15' being the least signficant. This tells us that a home having a waterfront view has the greatest impact on the price of the home while the square footage of the lot a home sits on has the least impact on the home price. The grade of the home, whether or not the home has a basement, and the number of floors are the next most impactful factors. It is interesting that bedrooms is so low on the list as one would usually think that homes with more bedrooms cost more. This is likely because the number of bedrooms can be arbitrary in a home. Two homes with the same square footage could have differing numbers of bedrooms and perhaps having larger but smaller rooms is more appealing to home buyers.
Running one final test on our model, using an OLS regression this time, we can check the p-values of each of our variables to make sure none of them are higher than 0.05 which would invalidate them being in our model
In [69]:
outcome = 'price'
predictors = df.drop(['price', 'price_log'], axis=1)
pred_sum = "+".join(predictors.columns)
formula = outcome + "~" + pred_sum
model = ols(formula=formula, data=df).fit()
model.summary()
Out[69]: